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I. ABSTRACT 

We present an overview of electronic device modeling using non-equilibrium Green function techniques. The basic 
approach developed in the early 1970s has become increasingly popular during the last 10 years. The rise in popularity 
was driven first by the experimental investigations of mesoscopic physics made possible by high quality semiconductor 
heterostructures grown by molecular beam epitaxy. The theory has continuously been adapted to address current 
systems of interest moving from the mesoscopic physics of the late 1980s to single electronics to molecular electronics 
to nanoscaled FETs. We give an overview of the varied applications. We provide a tutorial level derivation of the 
polar optical phonon self-energy Then, focusing on issues of a non-orthogonal basis used in molecular electronics 
calculations, we derive and the basic Green function expressions starting from their definitions in second quantized 
form in a non-orthogonal basis. We derive the equations of motion for the retarded Green function and the 
correlation function G^, and we derive the standard expressions for the electron density and the current that are in 
widespread use. We point out common approximations and open questions of which one finds little discussion in the 
literature. 



II. INTRODUCTION 



The applications of noncquilibrium Green functions have been extensive including quantum optics 4], quantum 
corrections to the Boltzmann transport equation Q , high field transport in bulk systems , and electron transport 
through nanoscaled systems. Our interest has been in this last category of electron transport through nanoscaled 
materials under a finite applied bias. Below, we review the work in this area and then provide tutorial derivations of 
the standard expressions. 

Over the last decade, non-equilibrium Green's function (NEGF) techniques have become widely used in corporate, 
engineering, government, and academic laboratories for modeling high-bias, quantum electron and hole transport in 
wide variety of materials and devices: III-V resonant tunnel diodes P, 0, 0i l l^ ITsL Q, 0, 0, ItR llsL ITol 

1^,11. 22, 23J, electron waveguides | 24 ]. superlattices used as quantum cascade lasers [25|, Si tunnel diodes l26l l27l . 
ultra-scaled Si MOSFETs [MlMlMEiLpl' Si nanopillars |3|Hj3|[, carbon nanotubes pi pi pi IbS [40, [4U 
^3, l4a. [43 , metal wires \4&L , organic molecules puL ISll |5^ 

JpJHl il ii El il 111 ifl il ll ^ 

Im. l65l 169. l67|. superconducting weak links f68], and raagnetic leads |38|. l69l l70l| . Physics that have been included 
are open-system boundaries (tJ, fuU-bandstructure 0, [iM 0, IH [t^I, bandtails [^l, the self-consistent Hartree 
potential 0, |^ [t^ , exchange-correlation potentials within a density functional approach [sl [sl [3| j 

acoustic, optical, intra-valley, inter-valley, and inter-band phonon^ scattering, alloy disorder and interface roughness 
scattering in Born type approximations 0, 0, Ull HI Hsl I20I Wl, '2E, "2^, 27], photon absorption and emission "2^, 



energy and heat transport Il6| . single electron charging and nonequihbrium Kondo systems j75t . l76li77i ..78. 79. 8Q>.81j. 
shot noise [11 E 113, A.C^ialsllsllsllsllsllsll. and transient response [Hil. Time-dependent calculations 
are described further in '9^ . General tutorials can be found in [^ . 

The general formalism for NEGF calculations of current in devices was first described in a series of papers in the 
early 1970s [^ 95] . The partitioning of an infinite system into left contact, device, and right contact, and 

the derivation of the open boundary self-energies for a tight-binding model was presented in |7l| . This theory was 
re-derived for a continuum representation in 93], tunneling through localized impurity states was treated in |94j . 
and a treatment of phonon assisted tunneling was derived in j95j . In 1976, the formalism was first applied to a 
multi-band model (2-bands) to investigate tunneling [9^ and diagonal disorder "9^, and in 1980 it was extended to 
model time-dependent potentials . 

The citation rate of the first article in the series fzij gives a good indication of the use of the NEGF formalism 
applied to electronic transport over the last 3 and a half decades. Of the 277 citations listed, 110 occur during the 
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years from 2000 to the present, 119 occur during the 1990s, 22 citations occur during the 1980s, and 26 citations occur 
during the 1970s. Over the last decade and a half, the citation rate has been steadily increasing. The motivation 
for the development of the NEGF tunneling formalism was the metal- insulator- metal tunneling experiments that 
received much attention during the 1960s 99]. The revival, or accelerated use, of the approach was motivated by 
the experimental investigations of mcsoscopic physics made possible by high quality semiconductor heterostructures 
grown by molecular beam epitaxy. In 1988, Kim and Arnold appear to be the first to apply the NEGF formalism 
to such a system, specifically, a resonant tunneling diode |^. As experimental methods progressed allowing finer 
manipulation of matter and probing into the nanoscale regime, the importance of quantum effects and tunneling 
continuously increased, and the theory adapted to address the current systems of interest moving from mesoscopics 
to single-electronics to nano-scaled FETs to molecular electronics. 

III. ONE DIMENSIONAL TRANSPORT THROUGH PLANAR SEMICONDUCTOR DEVICES 

A. NEMO 

Effort developing a quantum semiconductor device simulator for devices in which the potential varies along one 
dimension (ID) through planar semiconductors peaked during Texas Instruments' Nanotechnology Engineerin g pro - 
gram which ran from 1993-1997 and resulted in the the Nanoelectronic Engineering Modeling (NEMO) program |lO0j . 
The Nanoelectronic Engineering Modeling research program was initially conceived to model resonant tunnel diodes 
(RTDs) and, specifically, to discover the processes that determine the minimum valley current. The program included 
both theoretical development and software development. Furthermore, somewhat unique to this program, there was 
a strong experimental verification effort which consisted of the growth, fabrication, and measurement of hundreds of 
resonant tunnel diodes of varying geometries, layer thicknesses, and material systems forming a large test matrix. 

There were seve ral m ajor thrusts to the theory and modeling component: (a) a flexible treatment of the open system 
boundaries [ll ll^lToij . (b) charge self-consistency 0,123, Iz^ (c) incoherent scattering processes 0,123,1^123, 102J 
and (d) full bandstructure models . An example of a comparison between theory and experiment is shown in Fig. 
(Q. This is a simulation of a GaAs / AlAs RTD with a 5.66 nm well, 3.1 nm barriers, and 20 nm spacer layers at 
a temperature of 4.2K. The phonon peak is clearly visible. It rides on a background current resulting from interface 
roughness scattering. The scattering current is calculated self-consistently with the Hartree potential and a local 
density approximation (LDA) for the exchange-correlation. The magnitude of the scattering assisted valley current 
matches w ell th e experimental data. It was found that the Hartree potential significantly overestimated the intrinsic 
bistability 1 1 ( . Even with scattering and the LDA potential, the bistability is overestimated since it is observed in 
the simulations where it is not observed in the experiment. 

For room temperature RTDs, simulations indicated that the valley current was determined by thermionic emission 
through the second resonance for a number of RTDs in the test matrix. A comparison of full-band simulations with 
experimental data from the test matrix for Ino.47Gao.53As / AlAs RTDs and Ino.47Gao.53As / Ino.48Alo.52As RTDs 
is shown in Figs. Q and Q. The simulations model purely coherent transport and yet match well with the valley 
current of the experimental data indicating that at room temperature in these devices, the current is not limited by 
incoherent scattering. Purely coherent tunneling calculations we re ne ver able to match the valley current of 'notched 
weir RTDs with Ino.47Gao.53As / InAs / Ino.47Gao.53As wells |104| . These devices, which have the highest peak- 
to- valley current ratio of any semiconductor RTD, appear to have their room temperature valley current limited by 
incoherent scattering. 

As the NEMO program progressed, the simulation software was enhanced to model other devices and material 
systems to support the various experimental programs within TI and Raytheon. One of the important enhancements 
was the ability to model the Si / SiGe / Si02 material system. Fig. (0J shows close agreement between the experimental 
and simulated capacitance-voltage curve of a metal-oxide-Si (MOS) structure of Brar et al. [T"05j The approach taken 
for Si for quick design calculations as shown in Fig. |0J was to implement multiple de-coupled single-band models. 
In these models, Schrodinger's equation is solved independently for each band using Green function techniques. The 
quantum charge from each calculation is added, and the total is used in the Poisson solver. The quantum calculations 
and the Poisson calculation are iterated until convergence using a Newton-Raphson algorithm. The simulations of 
the C-V curve used a 4- independent band model consisting of one band for the 4 equivalent X- valleys, one band for 
the 2 equivalent X valleys, one band for the light holes and one band for the heavy holes. 

One of t he great successes of the NEMO software was its aid in the design of the first working Si / SiGe tunnel 
diode |106| |. There had been significant experimental effort building and testing devices with designs consisting of an 
8 nm Sio.5Geo.5 tunnel region sandwiched between Si. There was delta doping in the Si on cither side of the Sio.5Geo.5 
tunnel region and the heavy doping continued through the tunnel region with the n-p junction occuring in the center 
of the Sio.sGeo.s. Several iterations with NEMO lead to a new design with the length of the Sio.5Geo.5 tunnel region 



3 



cut in half from 8 nm to 4nm, and with the doping completely removed from the Sio.sGco.s tunnel region. The 
simulations showed that the delta doping on either side of the 4 nm tunnel region was sufficient to contain the electric 
field within the tunnel region and the band overlap between the n-conduction band and the p-valence band. The final 
design simulation is shown in Fig. |(SJ). A decoupled multiple single band model was used. In the strained Sio.sGeo.s 
region, the multiple bands are explicitly apparent since they are split in energy by the strain. 



Derivation of the Self-Energies 



Ever since its initial publication, there have been continued questions concerning the derivation of the self-energies 
described in the appendices of P|. For this reason, we will describe the details of the derivation of the polar-optical 
phonon self-energy given in Appendix A of 0|. 

We start with the general form of the electron-phonon potential, 



(1) 



where q lies within the first Brillouin zone, r is the electron coordinate, a and are the phonon annihilation and 
creation operators, and C/q contains the Fo urier transform of the electron-ion potential. Details of the derivation of 



this form of T4p can be found in Sec. 1.3 of [10 



To second quantize the electron coordinate r, we need to define the planar orbital basis and the electron field 
operators. The planar orbital basis for a zincblende (or diamond) lattice with 2 atoms per unit cell is 



(2) 



The anions sit on the Bravais lattice and the cations are offset by the vector v = |[111]. In the planar orbital basis, 
L is the layer index in the [001] direction where a layer includes a layer of anions and a layer of cations and the 
layer thickness is a/2, Rf is the coordinate in the x — y plane of the anion in layer L, Rf -I- Vt is the coordinate of 
the corresponding cation in the x — y plane, and k is a 2 dimensional wavevector in x and y. The indices and Ci 
label the atomic-like orbitals s, p^, d*^, s* on the anions and cations, respectively. With this basis, the electron field 
operators are then defined as 



V'(r)=E 



^ (r|aj, L, k)ca,,L,k + ^ {r\ci,L, k)cc,,L,k 



(3) 



where Cai,L.k is the destruction operator for state |ai,L,k) and Cci,L,k is the destruction operator for state |ci,L,k). 
The second quantized electron-phonon Hamiltonian is then. 



Hep - -y= ^ Uq (aq 



dVV'^(r)e^'5"'V'(r) 



k,k' L,L' 



^(a,L,k|e''ina',i',k')<^^c,,x^k' 

a, a' 

5](a,L,k|e'i'-|c',i',k')<i,kCc'.L'.k' 

a,c' 

-^(c,L,k|e^'i-|a',i',k')<i,kCa'.L',k' 



^ (c, L, kle^i nc', k')c^kCc',L',k' 



(4) 
(5) 
(6) 

(7) 
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To evaluate the matrix elements of e^'^ '" in Q - lO, we must expand out the planar orbitals in terms of the localized 
orbitals. The anion- anion matrix elements of Q become 

R,f ,R,J a, a' L.L' k.k' 

In the long wavelength approximation, the matrix element in |HJ is evaluated as 

(a,L,Rt|e^'i--|a',L',R;) = e''?^^^e''i*-^*<5L,L.5a,a'<5R,3J (9) 

where A = a/2 is the monolayer (anion plus cation layer) thickness. The phase e''' '' is evaluated at the position of the 
anion (Rt,LAz) and pulled outside of the integral. The orthogonality of the orbitals results in the Kronecker delta 
functions. Boykin has shown that, in general, r is diagonal in the empirical tight-binding basis, so that th e eva luation 
of the matrix element in is true in general and does not rely on the long wavelength approximation [108| . With 
the matrix evaluated as in ((SJ, term Q becomes 

a L k.k' Rt 



EEE^^'^'^^U.kCa.L.k-,, (10) 



a L k 



Evaluating the matrix elements as in Q, there will clearly be no matrix elements between anion and cation orbitals, 
so that terms ^ and © are zero. Following the same procedure to evaluate term Q gives 

EEE^"^^^"-'^''^^U,kCcx.k-,. (11) 

c L k 

The final general form for iJep is 



Hep - -j= ^ f/q (a, 



'q T- aLq 



L k 



/ . ^a,L,k^a,i,k-qt ^ ^ ^c,L,k^c,L,k-qt 



(12) 



To calculate the self-energy, we calculate the path ordered Green function expanding out the S-matrix in the 
interaction representation, 

GrL,.'L'(k;i,t') = -^{'Pe^^^^''^"'^^^c,,.LM{t)cl^^,,^{t')). (13) 

P is the path ordering operator, c is the Keldysh contour, and H' js) is the perturbing Hamiltonian, in this case, H^p. 
The brackets (• • • ) indicate the nonequilibrium ensemble average |lO9llll0j |. Expanding out the exponential to second 
order, gives two non-zero terms, the zero order term and the second order term. The first order term is zero since in 
the absence of interactions, (oq) = (aj^) = 0. 
The second order term is 

t)!^ E U^^U^^ E E E e-.^(-^+«-)e-.^(-^+«-) 

qi,q2 Li,L2 ki,k2 ctl,Q!2 

(■P ^ C^Sl ^ C^S2 (aqi(si) +a-qi(s2)) (aq2(s2) +aLq,(s2) 
''Qi.Li.ki (''l)''"! .ii.ki-qti ('*l)'^a2,i2.k2 (•*2)cQ2,L2,k2-qt2 {^2)Ca.L.Vi 

{t)^l',L'Mit')) (14) 

where is 1/2 if a is a cation orbital and zero otherwise. 

We now apply Wick's Theorem to factor the operators and evaluate them in the absence of interactions |lllL Ill2| . 
The phonon operators become 

(P (aqi(si) -|-aLq^(s2)) (aq2(s2) +aLq,(s2))) 

= (^aqi(si)a-q2(s2)) + (P'aq2(s2)aLqi(si)) OC (5qi,_q2. (15) 
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The Kronecker delta function arises because the operators are evaluated in the absence of interactions so that the 
phonon modes are conserved. 

There are 2 possible ways to factor the electron operators to give the rainbow diagram, Fig. H12(l . leading to the 
standard self-consistent Born approximation. Both factorizations are equivalent canceling the factor of 1/2 in Eq. 
(|14|l . There is also a bubble diagram identical in form to the Hartree self- energ y diagram. In bulk, due to momentum 
considerations, this diagram is exactly zero (p. 401 of Fetter and Walecka [H^. In a heterostructure, this term is not 
necessarily zero, and it has been investigated in detail by Hyldgaard et al. MTS]. We have ignored the bubble diagram 
and only consider the two equivalent factorizations leading to the self-consistent Born approximation. Choosing one 
of the two equivalent factorizations, the electron operators become 

= SkMl9aL,aiLi (k; t, Si )(5k2 ,ki -qn 5^1 Li ,02^2 C^l " Qtl ! > ■S2)'5k2-qt2,k5a2L2.a'L' 0^' ^2, (16) 

In (|16|) . a lower case g is used to indicate the bare Green function in the absence of the electron-phonon interaction. 
Since the operators are evaluated in the absence of phonon interactions, transverse momentum is conserved giving 
rise to the Kronecker delta functions. 

With the factorization of the electron and phonon operators, the second order term (|14|l becomes 



where 



q Li,L2 ai ,Q2 

dsijds2 [(Paq(si)a]j(s2)) + (Pa_q(s2)aLq(si)) 

9aL,a,L^(M^^ Si )5^^i^ (k - qt] Si , 33)5^2 L2 ^2, t') (17) 

^ ai — c, a2 = a 

— ^ ai — a, a2 = c (18) 
otherwise 

The general form of the path ordered Dyson's equation resulting from expanding Eq. 113|l is 

+ (ds (ds'Y, E 5rL,a,L,(k;i,s)E^^i^,„^i^(k;s,s')5r2L2,a'L'(k;s',t') (19) 

Comparing Eq. p9|l to Eq. (|17|) . we identify the self-energy as 

Sr,.„.2L2(k; si, s.) = 1 5: iC/ql^ e^^^^(^-^^+— ) 

q 

(Paq(si)a^(s2)) + (Pa_q(s2)aLq(si))] G^,i„„,i,(k - q*; si, S2) (20) 

In going from Eq. ((T7|) to Eq. (j^OJ), we have replaced the bare with the full that must be calculated self- 
consistently with the self-energy E^. This replacement of g^ with sums the whole class of rainb ow d iagrams an 
example of which is shown in the context of Hartree- Fock theory in Fig. 10.2 of Fetter and Walecka |ll2l |. 

The real-time self energies that we need, E^, and S-'^, are obtained from S-^ by considering the individual 
cases of si and S2 lying on different branches of the Keldysh contour. Or, one can simply use the relations of Langreth 
IQ.09J. For pedagogical reasons, we will evaluate S*^ by considering the time contour and not using the Langreth rules. 
S^(si,S2) is obtained from S^(si,S2) by considering the case of S2 > si. According to the path ordering operator 
P, the only way that S2 can always be greater than si is if S2 lies on the upper branch of the contour and si lies on 
the lower branch. Applying the path ordering operator, which places operators of later times to the left, to the time 
dependent terms of Eq. H20|l . we obtain 

(alj(s2)aq(si)) + (a_q(s2)a^q(si)) G< ,,^i^(k-qt;si,S2) (21) 



where 

h 



GaiLi,a2i2(^-q*;'5l'«2) = T (cl, ,i,_k_q, (S2 )Cai ,Li ,k-qt («! )) ■ (22) 
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The operators a and are evaluated in the absence of interactions. Therefore, their time dependence is 

aq(s) = aqe-'""^ 

The phonon ensemble averages are evaluated as 



(23) 



(24) 



where rtq is the average phonon number in mode q which, in equilibrium, is the Bose-Einstein factor. Using H24|l and 
(|23|l . the time dependent terms of (|22|l become 



^^-..,(..-..) ^ ^ 1) G< ^^,„^^^(k - q,; .1, S2) 



(25) 



In steady state, is a function of the time difference coordinate si — S2 which we Fourier transform into energy 
according to 



d{si — S2)e 



if (si-S2) 



+ (nq + 1) e-.(^i-^^) G< i^^„^i^(k - q,; Si - ^2) 



"qG< L,,„,L.(k - qt; S - fi^q) + (nq + 1) G< ^^_„^^^(k - q*; + fit^q) 



(26) 



Placing this back into Eq. (|20|l . we obtain the general form for the phonon scattering self-energy within the self- 
consistent Born approximation. 



A(Li-L2 + I'oi,a2) 



hqG< L„„,L2 (k - qt; - fi^q) + ("q + 1) G< ^^^^^^ - q, ; + tuj^)] 



(27) 



The evaluation of E-^ proceeds in exactly the same way with the difference being that si > S2 when evaluating the 
time-dependent terms of Eq. H2U|) . Therefore, the operator ordering is reversed in Eq. (|21|l . The result is 



, (k; £;) = 1 ^ |C/q|2 e.<?.A(L,-L2+.„„„2) 

q 

hqC^ Li,a2L2 (k - ^ + ft^q) + (^q + 1) G>^L„C..L,i^ ~ Clf, E ~ f^J^)] 

The relationship of the retarded self-energy T,^ to I]< and S> is 

E«(i, t') = e(t - i') [E> {t, t') - E< (t, t')] = -^e(^ - i')r(t, t') 

or, Fourier transforming the time difference coordinate to energy, 

X f dE' T(E') i , , 

Y.^{E) = P / — ^ T{E). 

^ ' 2tt E-E' 2 ^ ' 



(28) 



(29) 



(30) 



Since r(£') is Hermitian, the principal value integral is also Hermitian. This term gives rise to an energy renormal- 
ization but no relaxation or dephasing. Because the principal value integral is so difficult to perform numerically, this 
term is often ignored and E^(iJ) is approximated as ^V{E) where 



T{E) = ^[^>{E)-Y.<{E)\. 



(31) 



Another approach, the one implemented in NEMO and described in I*!, is to use the single-electron or low-density 
approximation. In this approximation S*" is set to zero since the inscattering into energy E is proportional to the 
electron density G< at energies E ± fiio. The electron density is assumed low and thus approximated as equal to 
zero. The outscattering E-* at energy E is proportional to the hole density G^ at energies EzL Hlo. Since the electron 
density is zero, the hole density is just the total density of states. Therefore G^ is replaced in Eq. I|28(l by —iA where 
A is the spectral function. All of the energy dependence of E^ is contained in the spectral function. Furthermore, 



from Eq. E= 



iF. Since G^ is related to A in exactly the same way as E^ is related to F, Eq. H3U|I converts 
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the spectral functions occurring in into retarded Green functions. The retarded self energy then has exactly the 
same form as Eq. (|28|l with replacing G^. This form of the self-energy does conserve current. 

With the general form of the self-energy 1)28(1 , we now consider the specific form of f/q corresponding to dispersionless 
polar optical phonons. For dispersionless polar optical phonons, o^q = tJo, fiq = riBifiujo) = \ j — 1), and 



e huj„ 



^oJ[q^ + qlY ^^^^ 

where Qo is the inverse screening length |ll4j . For numerical reasons, we prefer to have in Eq. H28|l be only a 
function of qt rather than k — qj. Therefore, we change dummy indices so that the expression for becomes 



(33) 



When it is appropriate to approximate G'^(qt; E) as being independent of the angle of qi in the x ~ y plane, then 
we can perform the angular integration analytically. This approximation would be appropriate near the conduction 
band edge of GaAs where the bands are spherical. The sum over q becomes the following integral in cylindrical 
coordinates, 



^G>{quE±h^,) 



27r 



£g-i9^A(Li-L2+i'ci.o 



de- 



|k-qP 



(|k-q|2-Kg2)^ 



(34) 



le 



The integral over 9 is 



2- fc2 

de 



k'^ + qf + ql- 2kqt cos 9 



[k^ + qj +ql + ql'- 2kqt cos 9) 
This integral is evaluated by making the substitution 



(35) 



(36) 



which converts the integral over 9 into a contour integral around the unit circle. The integral is then evaluated using 
the residue theorem. With substitution 1(36(1 . Eq. ((35(1 becomes 



kqt 



dz 



kqt 



bz + 1 [z 



bz 



(37) 



where 



b = 



fc2 + q2 + ql+ql 



kqt 



Evaluating l(37() with the residue theorem, we obtain 

1 



le{ql,q^,k'^'^ 



2-n- 



^{ql+qf + k^ + qlY - ^k^t M 
We write the integral over q^ in Eq. 1(34(1 as 

'"^ dq. 



qliql + qt + k^ + ql) 



■q^ + k^ + q^)^ 



cos [qzA{Li - L2+ J^ai.aa] hiql , ql , k^^ 



21 3/2 



2tt 



and the final form of the self energy becomes 



E<, , ,,{k;E) 



2tt 



1 



At 

47r2^ 



L- L 



,k\qt) 



"sGa,L;o',L'fe;^ + M + (nB + l)G>^^^,_^,(g,;£;-M 



(38) 



(39) 



(40) 



(41) 
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where 



L- L 



7 2 2 



tt/A 



dqz 



q,A{L-L 



tt/A cos 

dqz 



Vq^ + k'^ + qlY - Ak'^q^ 

qzA{L -L + u,)] ql{ql + + fc2 _^ ql) 



q\ 



F + g2)2 



(42) 



Note that the integrand in Eq. (|4HI is only a function of the magnitude qt = |qt|. The integral over the angle was 
done analytically. We keep the form of the full two dimensional integral because, in an effective mass model, we make 
the following change of variables. 



d^qt 
47r2 



/ deq^p2Di£qt) 

Jo 



dEzP2D{E-Ez 



(43) 



where Sq^ and are defined variables, E is the total energy, B is the bandwidth, p2D is the 2D density of states, 
and Emin is a minimum energy used in the numerical integration. Sq^ is the transverse energy using the dispersion 

relation in the flat band part of the injecting contact. For a parabolic dispersion, £q^ = where m\ is the effective 
mass in the injecting contact. Therefore, p2D — m* / {271^?) . The longitudinal energy is defined as E^ = E — e 



qt ■ 



B. Post NEMO ID NEGF Developments 

Since its delivery to the U.S. government in 1997, t he pu bhc availability of NEMO has been sporadic |ll5l| . Currently, 
NEMO is semi- publicly available with restrictions |ll6j| . Interesting post-NEMO developments of applications for 
NEGF to ID transport through planar semiconductor structures are the formulation of NEGF theory to model 
superlattices and quantum cascade laser s l25l and the application of NEGF theory to model the indirect phonon 
assisted tunneling in Si tunnel diodes |2^ l27|. 



1. Interband, Intervalley Tunneling in St 

The NEGF theory developed to model indirect phonon assisted intervalley, interband tunneling in Si tunnel diodes 
went back to the original theory of Caroli et al. described in the 4th article of the series . The original theory was 
rederived for localized orbital full-band models with the specific types of phonons important for interband tunneling 
in Si. A detaile d de scription of the theory is given in j27j . Calculations were performed both with a second nearest 
neighbor sp'^s* |ll7l | and a nearest neighbor sp^s*d^ planar orbital basis |ll8l |. The specific phonons of interest were 
the zone edge TA and TO phonons. From a theoretical and computational point of view, modeling interband tunneling 
in Si is interesting and challenging because it is primarily an indirect, phonon assisted process. Furthermore, it is an 
intervalley process in which an electron starts out in one of the 6 valleys near X along the A line and tunnels to the 
valence band at F. 

The band diagram of the device simulated along with the 6 valleys of the conduction band is shown in Fig. 
The tunnel current is dominated by the electrons in the 4 transverse valleys whose long axis is perpendicular to the 
direction of current flow. In the simplest picture, this is because the effective mass in the tunneling direction of the 4 
transverse valleys is 5 times lighter than that of the 2 longitudinal valleys, and their combined 2D density of states is 
5 times greater. There is a small coherent tunnel current that arises from electrons tunneling from the 2 longitudinal 
valleys of the conduction band directly to the valence band. This is possible since the transverse wavevector of the 
electrons in the 2 longitudinal valleys is centered at k ~ and can therefore be conserved in the coherent tunneling 
process. 

The c urrent modeled has 3 separate components. The coherent component is given by the usual expression P, l27l 

mini. 
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In Eq. (|44|l . lower case g symbols are used for the Green functions since they are the bare Green functions calculated 
in the absence of phonon scattering. The phonon assisted component of the current is given by |27| 



Eg. (|45() is Fermi's Golden Rule written in Green function form. It is a full-band version of Eq. (49) of Caroli et al. 
|95j | written for TO and TA zone edge phonons. To emphasize the intervalley process, subscripts and superscripts T 
and X are used in Eq. H45|) . The quantity a'^ indicates the left injected spectral function corresponding to Fig. 
Since it is injected from the left, it is in one of the X valleys of the conduction band. Similarly, indicates the right 
injected spectral function which is in the valence band near F. The integrations over the transverse momenta are 
performed around the respective valley minimum points. Eq. H45|l is evaluated individually for each type of phonon, 
TO and TA, and the results are added. The relative magnitude of the 3 components can be seen in Fig. [7| Note 
that the axis for the "Direct" current which is the coherent current is in mA whereas the axis for the phonon assisted 
current is in A. The magnitude of the coherent current is approximately 5 orders of magnitude smaller than the 
phonon assisted current. 

One issue that was theoretically investigated was to understand the effect of quantized states in the contacts. These 
formed in the delta-doped wells shown in Fig. The effect of these states was to give some structure to the current- 
voltage curves primarily in the negative differential resistance (NDR) region as is shown in Fig. [7| The inset of Fig. 
gives an intuitive explanation of why the structure occurs in the NDR region. The 2D channels are initially on. As 
they uncross, the channels get turned off giving structure to the I-V. Since the quantum wells are degenerately doped, 
inclusion of realistic broadening in the calculations completely wipes out any structure in the I-V )2gj . 

A comparison between the experimental and calculated TV curves is shown in Fig. |S1 The discrepancy in the 
magnitude of the experimental and calculated current (approximately a factor of 5) can be the result of several factors. 
The magnitude of the current depends quadratically on the magnitude of the deformation potential and exponentially 
on the width of the tunnel junction and the evanescant wave vector in the tunnel junction. A series of calculations 
were performed increasing the separation of the delta-doping planes and thus increasing the tunnel junction width. 
The magnitude of the peak current had a dependence of exp(— d/0.42) where d is the tunnel junction width in nm. 
For reference, the simplest analytic, parabolic, effective mass model results in a dependence of exp{—d/0.5). The 
dependence of exp(— d/0.42) agreed fairly well with the d epen dence of exp(— d/O.SS) extracted from reverse biased Si 
pn junctions from a later exhaustive experimental study |l2Clj |. 

The decay lengths are determined by the evanescent wavevectors in the bandgap. The full band calculations 
are compared with the parabolic effective mass model in Fig. |21 The full band calculations result in evanescent 
wavevectors in the gap that are smaller than that of the parabolic effective mass model. This is to be expected since 
the evanescent hole bands wrap around to an excited band above the conduction band and the evanescent electron 
bands near X wrap around to a hole band near X. The surprising thing is that in the simple analytical tunneling 
model, the smaller evanescent wavevector should result in a decay length larger than the effective mass model when 
in fact the current calculations show the opposite. This apparent contradiction is not understood. 



The most comprehensive application of NEGF theo ry to superlatt ices and quantum cascade lasers was presented in 
Lee and Wacker which built upon previous work |l2ll Il22l Il23j| . The difficulty with modeling quantum cascade 
lasers is that they are very long containing 40 or more periods of active and injector regions under large bias and a 
treatment of scattering is essential since the transport through the structure is all scattering assisted. Lee and Wacker 
reformulated NEGF theory to model these types of structures as shown in Fig. (TUl Their formulation of NEGF theory 
takes a radical departure from all of the other work discussed in this chapter ,2^] . All of the other work is based on 
a formulation of NEGF theory that has the work of Caroli et al. 0| as its foundation. In particular all of the other 
work partitions the infinite system into a "device" which is the system of interest and contacts which act as thermal 
reservoirs for electrons and the partitioning is done as first described in [tJI • In the formulation of NEGF theory by 
Lee and Wacker, there are no contact regions. Special "periodic" boundary conditions are used in which the Green 




+ 



[tr {a^(kr,S- file;)} f'^iE) (l - fiE - hw)) (ns(fi^) + 1) 
tr {alikr. E + hw)nB{hu;)] f'^iE) (l - f {E + hw)) 
tr {al{kr, E - hw)] (l - /^(i;)) f {E - nw)nB{nw) 
tr {a^(kr,S + hw)] (l - f'^iE)) f{E + hLu)inB{huj) + 1)]] . 



(45) 



2. Quantum Cascade Lasers 
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functions and self-energies satisfy 

G,^^yiE-qV) = G,,,iE) (46) 

In Eq. H46|l . i and j are the localized Wannier fmiction index within one period of the structure, E is the total energy, 
and V is the voltage drop across one period of the structure. Everything is periodic but shifted down in energy by 
the voltage drop across one period. 

Since there are no contacts in this theory, an expression for current density cannot be obtained by considering 
current flow from the contact to the device. Instead, an expression for the current density is obtained by evaluating 

The Hamiltonian is then broken up into an unperturbed part Hq and a scattering part H' . The unperturbed part 
gives rise to the standard current expression 

= / ^tr{t,,+iG< -t,+i,G<+J (48) 

where the matrices t and are in the basis of Wannier fuctions, the trace is over all Wannier functions and spin 
within a period, t is the off-diagonal Hamiltonian matrix coupling Wannier functions between periods i and i + 1, 
k is the transverse wavevector, and A is the cross-sectional area. We have used the fact that z is diagonal in the 
Wannier basis. This equation is equivalent in form to Eq. (18) and (22) of Caroli et al. l71| derived from the electron 
continuity equation, Eq. (20) of Lake et al. Eq. (20) of Svizhenko and Anantram '3l1, and Eq. H132|) derived 
later in this chapter. 

What is particularly interesting is that the scattering component of the Hamiltonian in Eq. H47|l results in a current 
term 

J_,, = _£. ^ 1 gtr {z [G<S^ + G«S< - S<G^ - S«G<] } (49) 

k 

This is the first time of which we are aware that such a term has been derived. In a basis in which z is diagonal, for 
example the Wannier basis or the empirical tight binding basis analyzed by Boykin 108], Eq. H49|l reduces to 

J«catt = ^E / ^tr{G<S-^ + G^E<-S<G-4-S«G<} (50) 

Under very general conditions, it can be shown that this term is equal to zero. See for example Sec. 2.4 of (6(1 or Eq. 
(67) of i. 



IV. TWO DIMENSIONAL TRANSPORT IN TRANSISTORS 



The primary objective for the development of quantum semiconductor device simulators in which the potential 
varies in two dimensions (2D) is the modeling of ultra-scaled Si based FETs |2^ |2^ [s^, |^ . As the device size and 
dimension increase, computational intensity increases, and the sophistication of the model must decrease. The 2D 
simulators treat the bandstructure in a multiple single-band model which includes 3 decoupled bands to account for 
the three non-equivalent valleys in the quantized channel for NMOS. The simulators have included Hartree quantum 
charge self-consistency, and open system boundary conditions at the source, drain, and gate. For full 2D quantum 
transport, empirical, elastic, incoherent scatteri ng mode ls have been implemented |33,|33|- The 2D problem has been 
further simplified by using mode-space models ^32, 'l24'| . 

True inelastic scattering has been included in a multiple subband model in which the subbands are only coupled 
through the phonon scattering |3lj |. Thus, the quantum calculations were reduced to multiple ID calculations greatly 
reducing the computational burden. This model was used to simulate the energy relaxation of electrons from dis- 
persionless deformation potential optical phonon scattering in Si G-type Intervalley scattering was included from 
phonons with energies of 12, 19, and 62 meV, The electron phonon interaction was treated in the self-consistent Born 
approximation, and the resulting quantum charge density was iterated self-consistently with a 2D Poisson's equation 
to convergence. 

Dispersionless, deformation potential optical phonon scattering results in a self-energy that is local in position. 
In ref. [Mj . the spatial scattering region was increased starting initially in the source, then moving through the 
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extension and channel, and finally into the drain to encompass the entire FET. The effect of inelastic scattering in 
various regions of the device is shown in Fig. Ijllll . The particular dual gate device modeled for this figure has a 10 nm 
channel, 1.5 nm channel thickness, and 1.5 nm oxide. The gate region lies between ±5 nm. As the scattering region is 
extended through the device from source to drain, the effect on the current is almost symmetric around the center of 
the channel. This indicates that for ultra-scaled FETs, scattering in the drain end is just as important as scattering 
in the source end. This contradicts the conclu sion of Lundstrom and Ren that "scattering in a short region near the 
beginning of the channel limits the on-current |l25j |." The reason for the importance of drain-end scattering in the 10 
nm FET is that the gate length is comparable to the scattering length so that hot electrons from the drain-end can 
be reflected back into the channel. As far as we are aware, this calculation was the first to predict these results. 

The largest scale simulations of a 2D FET have bee n perform ed by Jovanovic f30| . He performed full 2D simulations 
of both the MIT 90 nm and 25 nm bulk MOSFETs |30l Il26j . An elastic, empirical scattering model was used and 
calibrated to mobility data. For the 90 nm device, NEGF simulations, drift diffusion simulations, and experimental 
data were all essentially indistinguishable. However, the drift-diffusion and NEGF calculations diverge significantly 
for the 25 nm device. The differences arise from the threshold voltage shift due to the shape of the quantum charge 
distribution in the channel and tunneling through the gate region. 



3. Recursive Green Function Algorithm for G 



One of the significant algorithmic developments of the 2D codes was the use of the recursiv e Gre en function algorithm 
(RGFA) for ||2^ ■ The algorithm was first developed during the NEMO program |l27j and implemented in 
early prototypes of the NEMO code, but it did not make it into the final deliverable; only the recursive Green fuction 
algorithm for was used Q. In ref. [l^ the algorithm was used for coherent transport. In ref. [s^, the algorithm 
was used for transport in the presence of incoherent scattering. This was the first use of which we are aware of 
a recursive Green function algorithm used in the presence of incoherent scattering. It was developed to meet the 
intense computational demands of modeling the large 90nm and 25nm MIT bulk MOSFET devices. Below we give 
an overview of the algorithm. 

First, we note that for elastic incoherent scattering, can be factored into components injected from the left and 
right contacts. This can be seen from Eqs. (73) - (76) of ref. [l| or Eq. (13) of 30] . The empirical scattering model 
has a local, diagonal form for the self energies, 

S5 = '^r^,^/2 (51) 



(52) 



The form for S< conserves current. The general form for is 



G" 



G^S<G^ + G"E<^G 



(53) 



The first term on the right is from incoherent scattering, and the second term on the right is a source term due to 
injection from the contacts. The starting point for the RGFA for G^ is the equation for the left connected 



<A 



(54) 



The superscripts < and O are used to indicate whether the semi-infinite portion of g"^ , g^ , g^, or extends towards 
the left (source) or right (drain) side, respectively, erf j — ^^Gf^/ and GfJ = ifAf.j^ where the superscript s 
indicates source. Writing the equivalent equation for g'^^ and subtracting gives an equation for the left connected 
source spectral function 



^ i,i /IS I J- s . 



(55) 



and, similarly, the right connected source spectral function 



9L 



AS ,j. „>s J. 



(56) 
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The left connected source spectral function is connected to the source contact. The right connected source spectral 
function is not connected to either the source or the drain. It is only non-zero if there is incoherent scattering, i.e. if T 
is non-zero. With the right and left connected source spectral functions, the full source spectral function is calculated 
from 

Gh- (57) 
The above 3 equations must be iterated until convergence at which point current is conserved. 

V. THREE DIMENSIONAL TRANSPORT THROUGH CNTS, NANOWIRES, AND MOLECULES 

Three dimensional (3D) modeling has been applied to semiconductor nanowires, carbon nanotubes, and molecules. 
In all of the 3D applications of which we are aware, the transport is modeled as coherent throughout the nanowire, 
CNT, or molecule. This may be a good approximation in CNTs, but for electrons traveling through the HOMO or 
LUMO bands of molecules, the approximation is open to question. The 3D modeling can be roughly classified into 
two categories. In the first category, the Hamiltonian is constructed using an empirical tight binding model. For 
well- understood materials such as the common semiconductors, Si, Ge, and GaAs, and carbon nanotubes (CNTs), 
the empirical models can be hi ghly accurate since they are fitted to known properties of the materials such as band 
gaps and effective masses |ll8l I128ij . In fact, for semiconductors, the empirical models are more accurate than the 
ab-initio models. Carbon nanotubes have been extensively modeled using the TT-bond model [TtI [T2 i FTil ITp . In 
the TT-bond model, the basis consists of one pz orbital per atom perpendicular to the axis of the nanotube fl3^. The 
3D CNT problem has also been simplified by assuming cylindrical symmetry and using a mode-space model |il31i] . 
Applications of NEGF to 3D modeling of semiconductors have been more scarce. Both a full-band sp^d^s* model and 
a discretized single band effective mass model were used to calculate transmission coefficients through Si nanowires 
[33L I33. 135I I . No attempt has yet been made to combine the NEGF calculation with a self-consistent Poisson solver 
and calculate I-V characteristics of the Si nanowires under bias. 

In the second category, the Hamiltonian is constructed u sing ab initio models with den sity functional theory (DFT) 
and a localized orbital basis 111 SIM IS El El El El III 111 illMiSi^ category, 
NEGF has been integrated with existin g codes s uch as the quantum chemistry code Gaussian98 JfilL El 0, |3| , the 
computational materials codes SIESTA jia . Il33j and FIREBALL 134], and custom codes |37i i59l| . It is in this last, 
highly interdisciplinary category that we find some of the most interesting recent developments. 

The general field of molecular electronics has brought together quantum chemists, physicists, material scientists, 
and electrical engineers. The two primary fields from which theoretical approaches have been taken are solid state 
physics and computational quantum chemistry. Diagrammatic perturbation theory. Green's functions, and Dyson's 
equation are the foundati on of much of modern solid state physics theory, and it is the basis for non-equilibrium 
Green fun ction theo ry H llOTl Ill2| . Diagrammatic perturbation theory has had less usage in modern quantum 
chemistry |l35l Il36| . Conversely, non-orthogonal basis sets are often used in chemistry |l37l Il3^ , but are generally 
only cons i dered in physi c s by thos e working on computational material properties using a localized orbital basis 
|l39l ll4Ct Il4ll I142I I143I Il44l Il45j . One exception is Feynm an's derivation of the anticommutation properties of 
creation and annihilation operators in a non-orthogonal basis 1461. The derivations in the literature describing the 
integration of DFT with NEGF are extensive [37i lla. Isil 152. .54 .74 | , although there has been little discussion of the 
validity of DFT in non-equilibrium IH Il47l| . Recent work suggests that the breakdown of this approximation m ay b e 
the reason for the often observed large discrepancy between theoretical and experimental current magnitudes UA^. 
Also, the derivations make little use of s econd qu antized operators which are the natural language of Green function 
and diagramatic perturbation theory H Il07l Ill2j . 

Below, we provide tutorial level derivations of standard expressions starting from their basic definitions in second 
quantized form in a non-orthogonal basis for an effective single particle Hamiltonian. The focus will be on the new 
problems that arise when working in a non-orthogonal basis as compared to an orthogonal basis. We will point out 
open questions and standard approximations of which one finds little discussion in the literature. 



A. 



-lai-iA- 



+ t 



13 



VI. GREEN FUNCTIONS IN A NON-ORTHOGONAL BASIS 
A. The Non-Orthogonal Localized Orbital Basis 

We start with a non-orthogonal locahzed orbital basis {laj)} where the states correspond to the wavefunctions 
4>i{r) = {r\ai). We next define the biorthogonal basis such that 

(a,|/3,)=<5.,, (58) 

The overlap matrix is defined as Sij — {ai\aj). By inspection, the orbitals are written as linear combinations 

of the {|ai)} as 

1/3.) -EMt'^^'k. (59) 
p 

since this satisfies Eq. (|58ll . With the expansion (|59|l . we can calculate the inner product 

q . ' 

The last equality results from the fact that the inverse of a Hermitian matrix is also Hermitian. 
In our non-orthogonal basis {|ai)}, the identity operator is 

/ = E 1"^) ("^-1 = E = E (60) 

This is proven by sandwiching the identity operator between any two basis states, 

{ak\ai) = E,^2^l^['^"^]ij (^jWi) = [S]k,i 



We next define the creation, a|, and annihilation, a^, operators for the non-orthogonal basis states in the 

usual way such that aj creates a particle in the empty state \ai). In other words, starting from the vacumn state, |0), 
the single particle state \ai) is filled by applying the creation operator, 

al|0) = |a.). (61) 

The anti-commutation relations are 11461 



{a„a]}=5,,j. (62) 

and {tti, Oj} = |a|, = 0. We also define the creation, b"^ , and annihilation, b, operators for the bi-orthogonal basis 
states {lA)} such that 

bl\0) = lA). (63) 

Substituting the expansion of state \/3i) in terms of states \ctj), Eq. H59(l . into the right hand side of Eq. (|63|l . and 
then substituting in Eq. (|61|l for \aj), we obtain the transformations between the 2 sets of operators, 

p 

^3 = E[^"'],,p«f (64) 
p 

With these expansions and Eq. 162(1 , it is straightforward to calculate the anticommutation relations 

{b^A} = [S-']., (65) 
{a,,b]] = {at,6,}=(5,,, (66) 



14 



and {6„ b,} = {bib]} = {a,, b,} = {alb]} = 0. 

We also introduce the standard field operators, ip^{r) and V'(r), which are the creation and annihilation operators 
corresponding to the orthonormal position eigenstates |r). A particle in state |r) is created by 

^t(r)|0) = |r). (67) 

These have the standard anti-commutation relations, {ip^ {r'),tp(r)} = S{r — r'). By applying the identity operator, 
Eq. lEOl), to the right hand side of Eq. l|H7jl . 

|r) = ^K)[5-i]^^^.(a,|r) 

id 

= EN) 

hi 

= Y.\P,)^*ir) (68) 

j 

and then expanding the state \ai) using Eq. H61|l . we obtain the field operators in terms of the localized orbital 
operators, 

= (69) 
j 

Similarly 

V'(r)=^0,(r)6, (70) 



It is useful to introduce the 3 different kinds of matrix elements of any operator | 139| . We will use the identity 
operator as an example. The covariant matrix elements are 

= [S]^^^ EE (71) 

In other words, the overlap matrix is the covariant representation of the identity operator. The contravariant matrix 
elements are 

= [5-1] ^ (72) 

The inverse of the overlap matrix is the contravariant representation of the identity operator. The mixed representation 
of the identity operator is the usual identity matrix, 

ima,) = 6,., EE I) (73) 

where <5i,j is the usual kronecker delta function. These representations of the identity operator act like the metric 
tensor and serve to raise and lower indices, e.g., 

Y^h.P^'^^I^ (74) 

j 

From Eqs. H59|) and H64|l . we see that the states and the creation operators b] transform in the same way 

as contravariant quantities. If we were to use superscripts and subscripts to differentiate between the covariant and 
contravariant states, we could re-write Eq. (|59(l as 

p 

From Eqs. 169(1 and H70|l . we see that the field operators are contractions of a covariant and contravariant quantity 
and are thus independent of the localized orbital basis as they must be. We find this to be a useful check when 
performing derivations that any operator corresponding to an observable is a fully contracted product of covariant 
and contravariant quantities. 
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B. Non-equilibrium Green functions and correlation functions 



We will now derive expressions for the electron density and the current in terms of the relevant Green function 
quantities. Two correlation functions, 2 Green functions, and 1 spectral function are of particular interest. The 2 
correlation functions are defined by 



—I 



Gf,M = T{b]{t'Mt)) (76) 



GUtX) = -j-m)b]{t')). (77) 
The retarded and advanced Green functions are defined by 

G5(t,0 = Q(t-t')^{bmb]{t') + b]{t')h{t)) 

= e(i-i')[G>-(i,0-G< (78) 
G5(t,0 = e(i'-t)-^(6,(t)6](t')+fo](t')&.(i))- (79) 
The spectral function is defined by 

A.^At^t') = }^{b,{t)h]{t')+h]{t')b,{t)) (80) 

= z[G5(t,t')-G^,(i,0] (81) 
= z[G>.(t,0-G<.(i,0] (82) 

where the brackets (• • • ) indicate the non-equilibrium ensemble average and the creation and annihi- 

lation operators are in the Heisenberg representation. We will be concerned with steady state in which case the time 
dependence becomes only a function of {t — t') and is Fourier transformed to energy, e.g. 



G< - {E) = d{t- i')e*^(*-*'^/''G< - (t, t') (83) 



Note that these Green functions, being products of two contravariant quantities, are the contravariant representation 
of the Green function. It is this representation which is most useful for numerical computations. We will see that 
Gf^j{E) is given by the inverse of the matrix whose elements are the covariant representation of [E — H], i.e. the 
inverse of the matrix whose elements are {ai\E — H\aj) — ESij — Hi.j- This matrix is sparse and therefore desirable 
for numerical manipulation. 

The expression for the electron density is straightforward to obtain, and it is directly related to G^ . 

(n(r)) = (^t(r)v,(r)) 

= -*n^G< (t = <' = O)0*(r)0,(r) 

hi 

= -*E / (i?)^*(r)^,(r) (84) 

The expression for the current is less straightforward and we will defer its derivation until we have the equations of 
motion for G^ and . 

The equations of motion for the Green functions are derived by applying ihd/dt to expressions H76(l - (I78II and 
evaluating. 

^h:|-^(b]it'Mt))^{b]it')[b.^H]) (85) 
where the effective single particle Hamiltonian is 

^ = E^M^^&.- (86) 

i,3 
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The commutator in (|85|l is evaluated using the anticommutation relation Eq. H66(l . 



k,l 



k,l 



Placing this back into Eq. H85|) gives 



(87) 



(88) 



k.l 



Going through the same exercise for Gfj{t,t') results in the same equation of motion as for Gfj{t,t'). Applying 
ihd/dt to the retarded Green function, Eq. H78() . gives 

d 



hg-Gf,it,t') = 5it-t') [5-^],, +E 



k,l 



Multiplying through by S gives 



h J2 S^,lQ-Gf;^ {t, t') = S{t - t')6,^, + H^,iGf^j {t, t') 



or in matrix notation 



d 

ihS- - H 

ot 



G^{t,t') = lS{t-t') 



Fourier transforming to energy as in Eq. (|83|l results in 

[ES - H] G^(£;) = 1 

Similarly, Eq. H88|l becomes 

[ES-li]G<{E) =0 



(89) 
(90) 

(91) 

(92) 
(93) 



C. Boundary self-energies 

1. and 

Now one partitions the infinite system into a finite "device" region and "contact" regions [llllilllil. Tradi- 
tionally, one then includes the effect of the "contacts" on the "device" exactly using time dependent diagram mati c 
perturbation theory summed to all orders to obtain the Dyson's equations for each type of Green function [H FiIIm^ . 
Time dependent perturbation theory is for perturbations in the Hamiltonian which cause the states to evolve ac- 
cording to unitary transformations. We will see that now we need to also include a perturbation in the off diagonal 
elements of the overlap matrix. This, in effect, alters the basis states and thus the Hilbert space spanned by the 
states. Therefore, it is not clear to us how to incor por ate this into the formal none q uilibrium Green's function theory. 
What we will do is what many others have done |37l liM fsil Is^. fs^. 1^ Il39l | , which is to use matrix algebra to 
write an effe ctive Dyson equation for G^ to include the effect of the contacts. This is an application of inversion by 
partitioning |l49j| . 

We illustrate the concept with a concrete example of a periodic chain of atoms with nearest neighbor overlap of 
the atomic-orbital-like basis states. Atoms {— oo, . . . , 0} lie in the left contact. Atoms {1, . . . , N} lie in the "device" 
and atoms {{N -j- 1), . . . , oo} lie in the right contact. The matrix elements of the Hamiltonian group into intra-atomic 
subblocks Di^i and inter-atomic subblocks ti^i±i. The size of these matrices is equal to the number of orbitals per 
atom. In the matrix [£'S — H] , the matrix elements of the left contact couple to those of the device by the off-diagonal 
blocks 



ESio 



1-1,0 



(94) 



and 



ESo, 



(95) 



17 



The matrix elements of the right contact couple to those of the device by the off-diagonal blocks 

— tjv,Af+l = —tN.N+1 (96) 

and 

ESn+1M ^ ^N+l,N — -~iN+l,N- (97) 

For any atoms i,j G {1, . . . , N} in the device, we can write the standard Dyson equation for the retarded Green 
function using t in place of t. 

^1,3 — Sij + gi,ll^l,0'^Oj 

+ gfA'tAf,A' + lG^+ij (98) 

where is the "bare" Green function of the uncoupled device, and the energy argument is suppressed. We write a 
second Dyson equation for the exact G^^ and G^_|_j which cross the device-contact boundaries. 

Go"!, = go^to.iGf, (99) 

G^+lj = gw+l.Af + ltw+l.A'G^ J (100) 

where g^g is the "bare" Green function of the end block of the uncoupled left contact (also referred to as the 
"surface" Green function), and g^^_i is the "bare" Green function of the uncoupled right contact. These are 
generally calculated from the eigenmodes of the bulk contact regions EHHEliJ. Combining Eqs. (EHl), EHJ, and 
(|100|l gives the standard Dyson equation for the device Green function. 

~ gij + g«,l tl,Ogo, 0*0,1 ^i.j 

^1,1 

+ S^N^N,N+lgN+l,N+l''^N+l,N G^j (101) 

■^N,N 

The retarded self-energies resulting from coupling to the left and right contacts are, respectively, 

^lA = (ti,o - ^^Si,o) g^o (to,i - ^^So,i) (102) 

~ (tAr,Ar+l — ESn,N+i) Sn+1,N+1 (tw+l.AT — ESn+1,n) (103) 

We add the superscript B to indicate that these self-energies result from the open-system boundary conditions of the 
device, and they are not the result of dissipative processes within the device. Multiplying Eq. (|101|l by the N x N 
block matrix of the isolated device, [ESq — Hq] = [g^] , results in standard equation for the device Green function. 

[ESo - Ho - T,"^] = 1 (104) 



RB 
N,N- 



where has two nonzero blocks, Sf f and S 

At this point, we have derived the matrix equation for the contravariant retarded Green function, Eq. H92|) . using 
the Heisenberg equation of motion, Eq. H85|) , and then we used matrix algebra to determine the values of Gfj in the 
central device region in terms of the surface Green functions of the contacts, Eq. (|104|l . So far, we have not needed 
to know anything about time contours and nonequilibrium Green function theory. 

2. S< andG< 

To calculate physical observables such as the electron density and current, we need . To obtain G^, we need 
an expression for the self-energy E^^. This cannot be obtained without recourse to NEGF theory. However, as we 
noted above, perturbation of the basis states does not appear to be compatible with formal NEGF theory. We will 
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have to choose a form for S<-^ that appears to be consistent with the NEGF theory, that is physically reasonable, 
that satisfies the usual relations between Y.^, E^^, and S>^, and that reduces to the correct form in equilibrium. 

Once is fixed, Eqs. Hl()2|l and 1)103(1 . certain relationships between the various self energies must be satisfied 
such as 



i [S>^ - S<^] . (105) 

Furthermore, in equilibrium, is completely determined by from the relation 

S<'^ = ifT^ (106) 

where / is the Fermi factor. These two relations considerably limit the possible forms for S^^. One also finds the 
argument in the literature that since the contacts are in equilibrium, must equal i/lF^]^. However, S^f is a 

self-energy of the "device" consisting of matrix elements of atom 1 which is not in equilibrium. Therefore, we are not 
convinced of the validity of this argument even though it does give the correct answer. 

G"^ is obtained from the Dyson equation for which, in tu rn, is obtained from the Dyson equation for the contour 
ordered Green function by applying the rules of Langreth |7ll Il09| . The general form is 

G< =g<+g«I]«G<+g^S<G'^+g<S^G^ (107) 

For a single electron perturbation, H', the Dyson equation for the contour ordered Green function is |150) 



G^(s, ,s') = gP{s, s') + fds, fds2 g^(s, si)H'(5(si - S2)G^(s2, s') (108) 

J c J c 

where the integrals are along the Keldysh contour 0, Il09j . Since the perturbing potential H' is local in time, 
cx 6{si — S2), si and s? are always on the same side of the Keldysh contour resulting in the equation for and 
being, respectively | |7ll . il0a |. 

G^(i, t') = g«(i, t') +JdtiJ dt2 g"{t, h)n'5{h - i2)G«(t2, t') (109) 

and 

G<{t,t') = g<{t,t') 

+ Jdh J dt2s''{t,h)U'6{ti-t2)G<{t2,t') 

+ Jdh J dt2g<{t,h)Il'6{h-t2)G^{t2,t'). (110) 

In Eq. (|110|) . the term in Eq. (|107|l containing the self energy S< is missing because the local time dependence of the 
potential forces si and S2 in (|108|l to always reside on the same branch of the contour. We now want to check that 
our new energy-dependent effective potential t = (t — ES) is still local in time and that it does not introduce any 
memory effects which could give rise to a new term in the usual Dyson equation for G^. To check that t = (t — -ES), 
is still local in time, we need to consider Eqs. I|98l- llOOl) in the time domain. Writing down Eq. (|99|) in the time 
domain gives 



dti J dt2SoAt,ti 



d 

to4 + ihSo^i— ] 5{ti - 12) 



GUt2.t'). (Ill) 



Looking at the effective potential in the brackets, it is still local in time. Or, generalizing this to the Keldysh contour, 
ti and t2 could never be on opposite branches of the contour. Therefore, it appears reasonable that no new self-energy 
terms are introduced by the effective potential t and that by replacing t with t, G^ and can be derived from the 
usual Dyson equations for G^ as described in detail in ref. jlj]. 

Briefly, for completeness, we write down the Dyson equations for G^ including only the coupling to the left contact 
since the equations to include coupling to the right contact are identical in form and can be obtained by replacing 
subscript 1 with N and subscript with -I- 1. For any atoms z, j G {1, . . . , N} in the device, 

G< = g< + g.fiti,oG<^. + g,<iti.oGo^^. . (112) 
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The Dyson equations for the exact G^^ and G^^- which cross the device-contact boundary are 

= SofiiosGfj + g^Qto,iG^j. (113) 

and 

Gt = goVoaG^, (114) 
Substituting Eqs. pi4|l and 11131) back into Eq. (|112|) . gives the exact expression for G< of the device. 

G< = g< + g5 ti,og^oto.i G<^. + gfi ti,ogo<oto,i G^^,^- + g<i ti,ogo^oto,i (115) 



Multiplying Eq. pi5|) by the N x N block matrix of the isolated device, [ESq — Hq] = [g^] ^, and using [ESq — 
Holg*^ — 0, we obtain the equation of motion for the contravariant correlation function in the device domain 

[ESo - Ho - S^^] G< ^ ■E<^G^ (116) 
is obtained from H^^, Eqs. pU2|) and pU3|) . by replacing g^ with g<, i.e. 

^tf - (ti,o - ^Si,o) go<o (toa - ESo,i) (117) 

— i^N,N+l — ESn,n+i) Sn + i.n+i (tw+i,JV — ES^i+i^n) (118) 

The one question left is what is g^o(^) where for any 2 orbitals i,j associated with atom 0, jgi-^) ^ i{^\,fiia{^))'^ 
Previously, in an orthonormal basis, we say that creates the state \ajg) localized in the left contact which is in 
equilibrium by definition. Therefore, g^Q(£') is simply the spectral function times the Fermi factor. 



^ko{E)=^^{E~^iLW^{E) (119) 



where the spectral function a 



, and f{E — fi^) is the Fermi factor for an electrochemical potential 



of the left contact fiL- Now, however, creates the state 



n i 

which is extended. However, gQQ{E) is the noninteracting correlation function for which there is no overlap between 
the contact and device states, specifically, in our example, So,i = Si.o = 0. Therefore, is block diagonal with 
no matrix elements coupling the contacts to the device and the sum in (|12UII is nonzero only for atoms n in the left 
contact, n £ {— oo, . . . ,0}. Therefore, for the noninteracting Sqq, and big only create and annihilate states in the 
left contact which is in equilibrium so that 1119|) still holds. The equations for the correlation function and Green 
function are now formally identical to those derived previously for orthogonal basis sets 0. 

D. Current 

1. Standard Expression 

The derivation of the current operator generally starts with an expression for the local electron density which is 
then placed into the continuity equation to give an expression for the local divergence of the current [7l|. Instead 
of a local continuity equation, we will write a continuity equation for the total number of electrons contained in the 
orbitals that are elements of the "device." Note that we are defining the "device" as a set of orbitals rather than a 
region of space. The total electron number is 



N 

Nd = -if ■ 

a, 6=1 ia.jb 

-ihiT{SoG<{t,t)] (121) 



N 



20 



where a and b index the atoms in the "device," ia indexes an orbital of atom a, and the trace is over all states in the 
"device." We now take the time derivative of N]j as 



d_ 

di 



Nr 



-ih tr < Sn lim 



d_ 

di 



d_ 

dt' 



G<{t,t') 



(122) 



The equation of motion for G^(t, t') in the device is 

i?i|-SoG<(t,t')-HoG<(t,0 
at 



J dt, [S^^(t,ti)G<(ti,i') + S<^(t,ii)G^(ii,0] 



(123) 



and the conjugate equation is 



ih—G<{t,t')So-G<it,t')Uo 



J dti [G<(t,ti)S^^(ii,i') + G^(t,ti)£<^(ii,i')] 



(124) 



Subtracting Eq. H123|l from Eq. (|124|l . taking the limit t' t, and the tracing over all device states, we have 



dt 



tr{HoG<(t,t)-G<(i,<)Ho}- 

J dtitT{G<{t,ti)i:^^{ti,t) + G^{t,ti)i:<^{ti,t) 



-S«^(t,ii)G<(ti,t) 



(t,ti)G^(ii,i)}-0. 



(125) 



Eq. H125(l is the continuity equation for the total electron number within the orbitals that define the "device." In 
steady-state, — 0. The second term is also due to the cyclic invariance of the trace. The last term must also 
be zero, but we can break it up into contributions from the left and right contacts. These contributions will be equal 
in magnitude and opposite in sign. For our linear chain nearest neighbor example, the self-energies resulting from 
coupling of the device to the left contact are non-zero only for the orbitals of atom 1. Similarly, the self-energies 
resulting from coupling of the device to the right contact are only non-zero for the orbitals of atom N. 

Fourier transforming the last term of Eq. (|125|l , the particle current flowing out of the left contact into the device 

is 

/dJ^ 
— tr{G<i(i?)S^^f(i?) + Gfi(£;)S<f(i?) 

-S^^f (i?)G<i(i?) - i:<fiE)GUE)]} (126) 

where the trace is over the orbitals of atom 1. We regroup terms using the cyclic invariance of the trace to obtain 

^ = JSh'' ^^^^ (^^1 - ^^i) - - ^^'^^^ 

^ /|S*^-{rfi[/LAM+*G<,].} (127) 



This current expression was first written down by Meir and Wingreen |148| . In Eq. H127() . /l is the Fermi factor of 
the left contact, Tf^ — i (Sff — S^f ), and Ai i — i (G^^ — G^^) where, for example, Gf^^ is the sub-block matrix 

of G^ consisting of orbitals of atom 1. 

The integrand of Eq. (|127|l satisfies the following equality P, Il48j 

(128) 



r^i [/lAi_i + iGf_j] — to.iGJ^g — ii^o^Q i 
This is derived by inserting the Dyson equations for Gq and Gf q 

Go,i = go!oto,iGf 1 4- g^o^o.iGi^i 
Gfn = Gf iti_og(5;o + G<iti,og;^o 



(129) 
(130) 
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into the right hand side of Eq. H128() . Therefore, the general expression for the current Eq. (|127|l is also given by 

/dE 
— tr {to,iG<o~ti,oGo<i} (131) 

Since there is nothing special about atoms and 1 , we can say that Eq. (|127|l implies that the current flowing between 
any 2 atoms n and n + 1 in the chain is 

^ ^ j 27rfi*'^ {*",«+l^n+l,n - '^n+l,nGn,n + l } (132) 

or explicitly writing out the t terms, 



Eqs. (|131|l . (|132|l and (|133|l are identical to those obtained from empirical tight binding models |lll71|. 

The matrix elements (t — ES)/H can be thought of as a transition rate coupling orbitals from atoms n and n ± 1. 
It is not immediately clear that ES should enter into the dynamics. We will show below that these are the correct 
matrix elements which result in unitary transmission probabilities for a band of a periodic system. 

For coherent transport, G< = G-'^S^^G"^ and A = G^T^G'^ or writing out these expressions explicitly, 

Gti = Gf,, (z/Lrf Gl, + G^^^ (^/^rf G^,i (134) 

and 

Ai,i = G^^^r^jGf^^ + G^^jyTf jyG^ J. (135) 
Inserting the above two equations into Eq. (|127|l results in the usual form for the coherent current P, Ill9l Il48j . 

./ - y ^tr {ri,iG^^^r^,A,G^,i} ih - Jr) (136) 

This equation w as fir st derived by Caroli et al. [tJ] (see their Eq. (42)) and then later by Fisher and Lee |ll9j | and 
then may others |148| |. 

Once one has derived an expression for the transmission coefhcient, it is always wise to check that for a periodic, 
uniform system, the transmission is 1.0 independent of energy within a propogating band. We will now perform this 
check of Eq. H136() for a linear chain of atoms with one orbital per atom. The Hamiltonian matrix elements are 
{an\H\an) = e and {an\H\an±i) = t where n is the atom index. The overlaps are {an\an±i) = s and (q;„|q;„) = 1. 
Applying the Hamiltonian to the ansatz {i/j) = results in 

E^{t-Es)z-^+e+{t-Es)z (137) 

Substituting the phase factor z ~ e'*^ into l|137|) results in the E — k dispersion relation 

E - e + 2{t- Es) cos{k) = 0. (138) 

Since we are considering a uniform infinite chain of atoms, our device will consist of 1 atom, n — 1. The surface Green 
function = g^^ = g22 of both the left and right leads is given by the recursion relation 

gf =[E-e-{t-Es) gf [t - Es)] (139) 

Eq. (|139|l simply states that the surface Green function of the semi-infinite chain ending at site is the same as the 
surface Green function of the semi- infinite chain ending at site -1. Eq. (|137|l is a quadratic equation for z and Eq. 
(|139|l is a quadratic equation for gf . Writing out the solution for both equations, one finds that 

R ^ z ^ e'^ ^ cos(fc) + i sin(fc) 

t-Es t-Es t-Es ^ ' 

The boundary self-energy for both the left and right contact is 

E«s ^{t- Es) {t ~ Es) (141) 
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and 

= = = -2Ini E-^^ = -2{t- Es) sm{k) (142) 

The superscripts L and R indicate that the quantity is the result of couphng to the left or right contact, respectively. 
The exact Green function on site 1 is 

Gt - [E~e-2it-Es)gf {t-Es)]-' 

= [E - e-2{t- Es) (cos(/c) + i sin(fc))]~^ 
= [-2i{t- Es)sm{k)r^ 

" zi(ri + r«) = Ttb- 

The factor of 2 comes from the two identical self energies from the left and right leads. In the second to last 
line, the real part is zero from the dispersion relation Eq. H138f) . The transmission coefficient from Eq. H136f) is 
T = ^i^iG^iTi iGf^i which, for scalar quantities, is 

T = r^'|Gfi|' (144) 
where = = T^. Substituting Eq. (fl^ into Eq. (fH^ results in a unitary transmission probability. 

^^-^-1 (145) 

This demonstrates that the matrix elements t — £'S are the correct ones to use in the current equations. 

With these expressions for the current Eqs. (|127|l . H131() . H133|) . and (|136|l . we have re-derived NEGF theory for 
a non-orthogonal basis that is formally identical to that for an orthogonal basis. The equations for the correlation 
function G^, the Green function G^, and the current J, are now formally identical and can be obtained by simply 
replacing the off diagonal Hamiltonian matrix elements t with t — ES. All of the efficient numerical algorithms 
developed previously, such as the current expressions and recursive Green function algorithms, can be applied 0,[2l 



2. Direct Evaluation of Current Operator 

In deriving the standard current equation l|127(l . we started by calculating the current flowing out a specific set of 
orbitals rather than a region of space. Since our model has a specific basis, it would appear that we do not need to 
make this approximation, and we can directly calculate the surface integral of the current crossing a plane. However, 
we will demonstrate by example that for an incomplete basis such a straightforward calculation results in an expression 
for the current which violates particle conservation and does not reduce to unitary transmission for coherent transport 
within a band. 

For evaluating the surface integral of the current crossing a plane, we choose the plane between the atoms and 1, 
i.e. between the left contact and the device. The expectation value of the current crossing the plane is 

J = J ds- J(r) = ^ J ds- (i/;t(r)Vi/;(r) - [VV'^(r)] V(r)) (146) 

-^2 f 

= — ds - lim (V- V')G<(r,t;r',i) (147) 

2m J r'^r 

I ds. [0*(r2)V</).(r)-0,(r)V0;(r)] 
J2 G< (t, i) 1 ds • [</>;(r)V0.(r) - (/.,(r)V</>*(r)] (148) 



2m 



Eq. (|148|l is the one that we will use. We write it in the form of Eq. H147|) to show that for a complete basis, the 
value of the surface integral is independent of position along the wire. 
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We take the plane of integration to be the xy plane. Fourier transforming into the energy domain and taking d/dz 
of Eq. l(TT7|) results in 



d 



. J = J- [ ^ [ dx f dylim (-^^ ^] G<{v,v'-E) 



Using the equations of motion for (r, r'; E) within the device 



E 



2 m 



V{v) 



G<{y.y'-E) = 



(149) 



(150) 



Eq. (|149|l becomes 



E+- V{r' 

2m 



dE r . r . ( 



G<{r,r';E) = 



dz 2m J 2Tih 
Integrating by parts using 



dx I dy lim 



52 



92 



92 



r'^y \dx'^ dy'^ dx'^ dy'' 



G<{v,v'-E) 



|G<(r,r;i.).lhn^(^ + |,)G<(r,r';i.) 

Eq. (|152|l is equal to zero. This is perhaps easier to see by re- writing Eq. (|152|l back in the form of Eq. H146() . 
We now return to Eq. p48|) which we write out explicitly as 



(151) 



(152) 



(153) 



J = Y.^U^^t) Jdxjdy 



2m 



d 



d 



c^*{v) — Uv)-Uv)—J*{v) 



dz^' 



(154) 



We evaluate (| 154(1 at z between sites and 1. The only nonzero contributions come from the orbitals that cut the 
surface between atoms and 1. In our nearest neighbor example, only the orbitals of atoms and 1 contribute to 
the integral. The integral in (|154|l contains 4 terms, 2 intra-atomic terms and 2 inter-atomic terms. We will first 
consider the inter-atomic terms since these 2 terms will result in an expression formally identical to Eq. H131|l but 
with different matrix elements. 
The inter-atomic terms are 



(155) 



where is the matrix element ^ dx ^ dy {^-^^ (j)*{r)-^(j}j{r). 
We define the current matrix element as 



and Eq. (|155|l becomes 



= tr {Jo,iG<o(t,t)+Ji,oG<i(i,i)} 

= / {■^OAG<,iE)~jl,G<,{E)} 



(156) 



(157) 



where in the last line we used jj = — Ji.o. 

Eq. p57|l is a general expression for the current. We can move the plane to cut between any two atoms of the 
chain, so that we have an expression for calculating the current at any point in the device. For example, cutting 
between atoms n and ?i + 1 , the current would be 



(158) 
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From now on, we will work in the energy domain and save on notation by defining Jt{E) as Jt{E) = 

tr |jo,iG^Q(i5) — jj ^G^;^(-E)|. We will also use the matrix notation of the last hne of Eq. (|157ll . The bold 

quantities indicate matrices the size of the localized orbital basis per atom and the subscripts indicate from which 
atom the orbitals come. The trace is over the orbitals equivalent to the in Eq. (|157|l . 

To bring Eq. p57|l into a form comparable to Eq. H127|l or Eq. (| Infill , we begin by substituting the Dyson Equations 
for G<i and G<o Eqs. ifT^ and lfnO|) into Eq. l(T77|) to obtain 



Jt{E) — tr |jo,iGf ^ti^ogao + JoaGJ^^jti^ogo^o ^ Jo.igo^o^oaGJ^^ — jj ^^gJ-QtciG^^j 



(159) 



Eq. (|159|l is the equivalent of Eq. I|127|) . If the Js were replaced by ts, then Eq. (|159|l would be identical to Eq. 
(|127|l . Grouping the terms of Eq. I|159|) with G^j^ together and using the cyclic invariance of the trace, Eq. H159|) 
becomes 



Jt{E) = tv li 



(jo,igo!oto,i) - (jJ,igo^oto,i) 

+tr |ti,ogiJ;oJo,iG!-f 1 — Jaigoio^i^oGf ^| 
tr 7i^^i iG<i + tr {ti,ogo<o Jo.iGf i " AA'^i-^^ti ] 



where we used tj^ = ti_o, [go^ol^ ^ gc^O' ^'^'^ defined 



7i,i = i 



Jo.igaotoa) - (Joigaoto.i 



(160) 



(161) 



Eq. (|160|l is exact, valid with or without incoherent processes in the device. 
The 2 intra-atomic terms from Eq. (|154|l result in a current term 



— tr {Jo,oGo<o(^^)+Ji4Gm(^^)} 



(162) 



In equilibrium when G^ = i/A, the intra-atomic terms are individually since for our real orbitals, A(t, t) = S ^ 
is symmetric and J is antisymmetric. Gq Q(i?) is the exact G*^ in the left contact, which is, since the contact is by 
definition in equilibrium, G'^q{E) = iJlKq^qIE) |l5lj| . Therefore, the first term of Eq. p62|) is always zero for a basis 
of real orbitals. 

For coherent transport, we can br ing t he sum of Eqs. (|16U|I and Eq. (|162|l into the Fisher-Lee form for the 
transmission coefficient. Substituting |151| 



^1,1 — ^IlGi^C^iG^i + ijRG^j^T^j^Gfj-^ 

into Eq. p6U|l gives 

Jt{E) = -tr{7fiGf;jvr^jvG;^.i}/i? 



-|-tr |— 7'/^jG^^r^jG^^ + i ti^oao,oJo,iG{^^ — jj ;^ao,oto,iG^j | 



(163) 



(164) 



where in the second line we used the relation g^Q = i/^ao.o. Substituting Eq. p63|l into the intra-atomic current 
term Eq. ifTC^ gives |l5lj 



(165) 



The total coherent current is the sum of the intra-atomic Eq. p65|) and inter- atomic Eq. (|164|l contributions. 

In equilibrium, when Jl — fn, the current must be 0. Furthermore, we can, in principle change /l and fu 
independently of the transmission coefficient. Therefore, the term proportional to /l must be equal in magnitude 
and opposite in sign to the term proportional to /r. Both terms are alternative expressions for the transmission 
coefficient. The term proportional to fa is in the standard Fisher-Lee form so that we can write 



— tr {{ri, - zJi.i) G^^^rf ,^G^,i} ih - Ir) 



(166) 
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We note that the other term of Eq. (| 164(1 is the more computationally efficient form equivalent to Eq. (53) of ref. jj]- 
Also, I is Hermitian and the product Gf n^n i Hermitian. Therefore, the trace tr {7i.iGf n^n i} 
is real. Conversely, Ji^i is antisymmetric so that i tr jJi^iGf^F^ is also real. Eq. (|165|l is identical to the 

expression that is universally used l^l S El El El Ei En. fH^ with the replacement of Fi i with {"Ji i — iJi.i). 

With current expression we investigate the consequence of an incomplete basis by evaluating the value of the 

transmission probability for a periodic, uniform system as we did for the "standard" current expression Eq. (|136|l . 
Using the same single-band model described immediately preceding Eq. p37|l . the intra-atomic term Ji_i is 0. The 
factor is 

7-^ = -2Im JQ,igf{t - Es) = -2 Jo,i sin(A:). (167) 



The transmission coefficient Eq. (|166|) is then 



T=^. (168) 



The transmission probability of Eq. (|168|l can only equal unity if Jq.i = t — Es. Since t — Es is energy dependent, 
this is clearly not possible. Therefore, this "exact" evaluation of the current leads to qualitatively incorrect results for 
an incomplete basis. With an incomplete basis, one should derive the current expression from conservation laws that 
are consistent with the Hamiltonian and basis such as Eqs. H121|) - (|125f) . We note the strong parallel here with the 
discussions in the literature concerning t he form of the mom entum matrix clement in empirical tight binding theory 
in which the explicit basis is not known jlOSl Il52l Il53l Il54j . 



3. Equivalence of the Standard and Exact Current Expressions in a Complete Basis 

The argument for the equivalence of the two current expressions, Eqs. H127|l or H131|) and Eqs. I|157|) plus H162|) is 
base d on particle conservation by writing the current as the time derivative of the electron number in the left contact 
|l55j . For lack of better expressions we will refer to the current in the "standard" expressions Eqs. ((127|l or 1)131(1 
as the "tight-binding current" {Jtt)- We will refer to the current derived from the direct evaluation of the current 
operator Eq. 1(1471) as the "real space current" {Jz)- 

For the real space current, Eq. 1(147(1 is derived from the continuity equation, and it is equal to the negative time 
derivative of the electron number in the left contact (for the surface chosen between atoms and 1). It is given by 

J. ^ ~^ = -g-J dx J dy dz (^t(r)VXr)) 

= J'^'^J'^yJ'' dzcb*{r)Mr) (169) 

where is the electron number of the left contact and zq lies between atoms and 1. 

The tight binding current Eqs. 1(127(1 or 1(1311) were derived by considering the time derivative of the electron number 
of the orbitals in the "device." It can also be derived by considering the time derivative of the electron number of the 
orbitals in the left lead. Although this seems physically obvious, it is not straightforward to show, therefore below, 
we write down the main steps. 

The electron number in the left contact is given by the sum over all of the orbitals of the left contact. 




= -ihti{SoG<{t,t)} (170) 

where the trace is over all orbitals of the left contact. We need an equation for the exact of the left contact in 
terms of the exact Green functions of the device. Mimicing the steps used to derive of the device Eqs. ((112(1 - 
((116(1 we write the Dyson equations for the left contact. For any atoms i,j G {— oo, ... ,0} in the left contact, 

^tj = + gfotoaGfj + gj^otoaG^^j. (171) 
The Dyson equations for the exact G^^ and G^^ which cross the device-contact boundary are 

G<^. = G{^iti,ogo1, + G^iti^og^^,. (172) 
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and 

Gt = GliKog^.j (173) 
Substituting Eqs. H173|l and l|172(l back into Eq. H171(l . gives the exact expression for G< of the left lead. 

G< = g< + gfotoaG^^iti^ogo^, + g5to,iG<iti,,g^^. + g<oto.iG^;iti,og^, (174) 
Muhiplying ((TTljl on the left with g^"^ = [ESa - Ho] results in 



[SSo — Ho]j G^,^- — 5i,o to,iGf iti^o g^j + 5i^o to,iGf jti^o go^j- 



•^0,0 



(175) 



In the time domain we have 



zh^SoG<{t,t')-IloG<{t,t') = J dh [S«^(t,ii)g<(ti,0 + S<^(t,ii)g-^(ti,t')] 



and 



(176) 



-ih-^G<it,t')So-G<it,t')Uo= J dii [g<(t,ti)S^^(ii,t')+g''(i,ii)S<^(ti,t')] (177) 
Subtracting Eq. p76|l from Eq. (|177|1 . taking the limit t' t, and the tracing over all contact states, we have 

^+tr{HoG<(t,i)-G<(t,i)Ho}- 



J dhtT{g< {t, tl)S^^(tl, i) + g^(t, il)S<^(tl, i) 

-s«^(t,ii)g<(ti,i) - i:<^{t,h)g^{h,t)} = 0. 



(178) 



The second term in Eq. H178|l is zero due to the cyclic invariance of the trace. Writing out the last term we get that 
the current flowing out of the left contact is 



Jtb — 



tb 



dt 



dE 



- I ^;^tr {to,iGi iti,ogo,o + to,iGi,iti^ogo,o 



(179) 



Cyclically permuting the terms under the trace, Eq. (|179|l is identical to Eq. (|126|l derived by considering the time 
derivative of the electron number of the device. 

Now that we have shown that both currents are equal to the negative time derivative of the total electron number 
of the left contact, we take the long time average of the current. 



{Jz 



1 
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dtJ, 



[Ni{r)-N^{-T)] /27 



(180) 



Similarly 



{Jtb) = ^ f_ dtJtb = - K(r) - N,^,i-r)] /2r 



(181) 



The only difference between these two expressions lies in the orbitals near the plane between atoms and 1. The 
integral defining contains part of the orbitals of atom and part of the orbitals of atom 1. The integral defining 
contains all of the orbitals of atom and none of the orbitals of atom 1. We can say that the difference between 
(Jtb) a-nd ( J^) must be bounded by 



a,b—0 i,j 



/2r 



(182) 
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where a and b label the atom index and i and j are the orbital index. If the quantity within the vertical brackets 
is finite, then the long time average of the two current expressions must be identical. We note that a rigorous 
m athem atical proof requires careful construction of the scattering states an example of which is given in Appendix A 
of 1156*1 ■ and also care in taking the appropriate long time limits for both the scattering state and the adiabatic turn 
on of the device-contact coupling. 

To conclude this section, we re-state the fact that the "standard" current expressions are the correct ones since they 
result in a unitary transmission probability within the band of a periodic system. The "standard" current expressions 
Eqs. H127|l and (|132|l are identical in form to the current expressions derived for an orthonormal basis. The only 
difference is that the Hamiltonian matrix element t is replaced everywhere by the effective Hamiltonian matrix element 

Even with an explicit basis, one should derive the current expression from conservation laws that are consistent 
with the Hamiltonian and basis such as Eqs. (|121|l - (|125|l . A direct evaluation of the real space current operator in 
an incomplete basis will lead to a breakdown of the particle conservation law and qualitatively incorrect results. For 
a complete basis, the 2 approaches should give the same result. 

VII. CONCLUSION 

In summary, the use of the NEGF approach to quantum electron transport developed in the early 1970s is continuing 
to increase. So much so that one could argue that it is the most heavily used approach for modeling quantum electron 
transport through nanostructures ranging from semiconductors to carbon naontubes to molecules. 

Currently, the heaviest usage is in the area of molecular electronics in which the NEGF approach is used to model 
current-voltage characteristics of molecules. All of the NEGF applications in this area, of which we are aware, 
model coherent transport through the molecule. The majority of the reported implementations combine DFT with 
NEGF often integrating the NEGF algorithm with existing ab-initio commercial quantum chemistry or computational 
materials codes. These codes use an explicit non-orthogonal localized orbital basis. The formulation of the standard 
NEGF approach to open system boundaries in a non-orthogonal basis raises questions that have not been addressed 
in the literature. 

To address those questions, we have re-derived the standard NEGF theory in a non-orthogonal basis using the 
second quantized formalism that underlies the theory. This has allowed us to explore some of the approximations 
that are commonly made, but of which, one finds little discussion. The most fundamental approximation lies in the 
compatibility of NEGF based on adiabatic perturbation theory and unitary evolution of states with changes in the 
basis states and corresponding Hilbert space. This issue does not arise if one is only interested in since NEGF 
theory is not required to get G^. It can be obtained from the Heisenberg equation of motion and matrix algebra 
as shown by the derivation of Eq. (|HJ4|I . In nonequilibrium, one must obtain from the contour ordered Dyson 
equation, which, in principle, does not appear to be compatible with perturbations in the basis states. However, once 
Yi^^ is fixed by the derivation of G^, the options for and thus G^ become limited, and the standard expression 
for S*~^ Eqs. H117|) and H118|) found in the literature appears reasonable. 

The second approximation underlies the derivation of the "standard" current expression Eq. H127I136() . We were 
unable to find a form for the local electron density from which one could derive these expressions from a local continuity 
equation. Instead, we wrote a continuity equation for the total number of electrons in a set of orbitals that define the 
"device" or left contact rather than in a specific region of space. In the coherent limit, the resulting expression for 
the transmission probability gives the correct unitary value for a periodic system. 

A direct evaluation of the real space current operator in an incomplete basis will lead to a breakdown of the particle 
conservation law and qualitatively incorrect results. For a complete basis, the standard approach and the direct 
evaluation of the real space current operator should give the same result. 
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Figure Captions 

1. Example calculation including incoherent scattering from acoustic phonons, polar optical phonons and inter- 
face roughness. The simulations are compared against experimental data. Reprinted from 21] . Copyright 1996 IEEE. 

2. Test matrix of Ino.47Gao.53As / AlAs RTDs. The collector barrier is increased by one monolayer for the 3 RTDs 
from left to right. Forward and reverse bias current voltage curves are shown and overlaid on experimental data. 
Reprinted from Copyright 1997 IEEE. 

3. Test matrix of Ino.47Gao.53As / Ino.48Alo.52As RTDs. (a-c) The well thickness is increased, (c-f) the undoped 
spacer layer is increased. Reprinted from 22]. Copyright 1997 IEEE. 

4. Experimenta l dat a and simulation of the capacitance - voltage curve of a MOS structure with a 3.1 nm oxide. 
Reprinted from Copyright 2001 IEEE. 

5. Simulation of band profile of the Si / Sio.5Geo.5 tunnel diode. Reprinted from |l06j ]. 

6. Band profile of the delta-doped Si tunnel diode biased at O.IV. Inset, 6 X valleys of the conduction band. 
Reprinted from p^ . 

7. Three current components: TA phonon assisted, TO phonon assisted, and direct or coherent current. Inset, 2D 
electron and hole dispersions. Reprinted from p^ . 

8. Experimental I-V from ref. |157| corrected for 12 fl series resistance overlaid on calculated I-Vs with different 
energy broadening in the contacts as shown in the legend. Reprinted from |23] . 

9. Comparison of the real and imaginary dispersion relations calculated (a) from the parabolic single band model 
using the transverse conduction band mass of 0.19 mo and from the full-band model in the transverse mass direction 
of the A4 conduction band valley and (b) from the single band model using the light hole mass of 0.16 mo and from 
the full-band model in the valence band in the (001) direction. The horizontal axis to the left of is imaginary k and 
to the right of is real k. Reprinted from pTi ]. 

10. Quantum cascade laser structure modeled in Reprinted with permission from Copyright 2002 American 
Physical Society. 

11. Drain current and channel barrier height versus extent of scattering region (starting from source at -20 nm). 
Inset: I^i vs. VdS curves for Vq — 0.6 V. Reprinted with permission from 31]. Copyright 2003 IEEE. 

12. Self-energy diagram for self-consistent Born approximation. 



13. Keldysh time contour. 
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Experiment 
Forward Scatter 




FIG. 1: Example calculation including incoherent scattering from acoustic phonons, polar optical phonons and interface rough- 
ness. The simulations are compared against experimental data. Reprinted from |21| . Copyright 1996 IEEE. 




FIG. 2: Test matrix of Ino.47Gao.53 As / AlAs RTDs. The collector barrier is increased by one monolayer for the 3 RTDs from 
left to right. Forward and reverse bias current voltage curves are shown and overlaid on experimental data. Reprinted from 
[23. Copyright 1997 IEEE. 
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FIG. 3; Test matrix of Ino.47Gao.53 As / Ino.48Alo.52As RTDs. (a-c) The well thickness is increased, (c-f) the undoped spacer 
layer is increased. Reprinted from j22l] . Copyright 1997 IEEE. 
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FIG. 4: E xperimental data and simulation of the capacitance - voltage curve of a MOS structure with a 3.1 nm oxide. Reprinted 
from lll^. Copyright 2001 IEEE. 
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FIG. 6: Band profile of tlie delta-doped Si tunnel diode biased at O.IV. Inset, 6 X valleys of the conduction band. Reprinted 
from |2^. 
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FIG. 7: Three current components: TAphonon assisted, TO phonon assisted, and direct or coherent current. Inset, 2D electron 
and hole dispersions. Reprinted from |2q|. 




FIG. 8: Experimental I-V from ref. |157| | corrected for 12 n series resistance overlaid on calculated I-Vs with different energy 
broadening in the contacts as shown in the legend. Reprinted from |2^ . 
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FIG. 9: Comparison of the real and imaginary dispersion relations calculated (a) from the parabolic single band model using 
the transverse conduction band mass of 0.19 mo and from the full-band model in the transverse mass direction of the X4, 
conduction band valley and (b) from the single band model using the light hole mass of 0.16 mo and from the full-band model 
in the valence band in the (001) direction. The horizontal axis to the left of is imaginary k and to the right of is real k. 
Reprinted from |27|. 
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FIG. 11: Drain current and channel barrier height versus extent of scattering region (starting from source at -20 nm). Inset: 
Id vs. YdS curves for Vg = 0.6 V. Reprinted with permission from Isi!]. Copyright 2003 IEEE. 
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FIG. 12: Self-energy diagram for self-consistent Born approximation. 



FIG. 13: Keldysh time contour. 



